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Abstract 

We summarize the theoretical approach to the solution of the NNLO DGLAP equations 
using methods based on the logarithmic expansions in x-space and their implementation 
into the C program Candia 1.C0. We present the various options implemented in the 
program and discuss the different solutions. The user can choose the order of the evolution, 
the type of the solution, which can be either exact or truncated, and the evolution either 
with a fixed or a varying flavor number, implemented in the varying-flavour-number scheme 
(VFNS). The renormalization and factorization scale dependencies are treated separately. 
In the non-singlet sector the program implements an exact NNLO solution. 
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PROGRAM SUMMARY 



Program Title: Candia 

Journal Reference: 

Catalogue identifier: 

Licensing provisions: none 

Programming language: C and Fortran 

Computer: all 

Operating system: Linux 

RAM: In the given examples, it ranges from 4 MB to 490 MB 

Keywords: DGLAP evolution equation, parton distribution functions, x-space solutions, QCD 
PACS: ll.10.Hi 

Classification: 11.5 Quantum Chromodynamics, Lattice Gauge Theory, 11.1 General, High Energy 
Physics and Computing. 

Nature of problem: 

The program provided here solves the DGLAP evolution equations for the parton distribution func- 
tions up to NNLO. 
Solution method: 

The algorithm implemented is based on the theory of the logarithmic expansions in Bjorken-x space 
Additional comments: 

In order to be sure to get the last version of the program, we suggest to download the code from our 
official website of Candia is http://www.le.infn.it/candia. 
Running time: 

In the given examples, it ranges from 1 to 40 minutes. The jobs have been executed on an Intel Core 
2 Duo T7250 CPU at 2 GHz with a 64 bit Linux kernel. The test run script included in the package 
contains 5 sample runs and may take a number of hours to process, depending on the speed of your 
processor and the amount your RAM. 

LONG WRITE UP 
1 Introduction 

Perturbative predictions in QCD are going to be essential for the discovery of new physics at 
the LHC, the new hadron collider at CERN, and for this reason the determination of the QCD 
background for important processes requires the analysis of cross sections at higher perturba- 
tive orders in an expansion in the strong coupling constant a s . By now, the level of theoretical 
accuracy reached in the calculation of several hadronic observables for LHC studies is rather 
impressive, and this has been possible thanks to the development of new perturbative tech- 
niques which have allowed to move from previous next-to-leading-order (NLO) determinations 
of several key processes, to the next order in accuracy, which is the next-to-next-to-leading 
order (NNLO). These computations of the hard scatterings need to be accompanied by the 
corresponding NNLO DGLAP evolution in order to be phenomenologically and theoretically 
consistent. Some available codes that deal with the DGLAP evolution are Pegasus [1], based 
on the use of Mellin moments, Qcdnum [2] and Hoppet [3], both based on a discretization of 
x and fj/p. 

There are several issues that need to be addressed when we move to this perturbative 
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order. One of them concerns the size of the corrections, which are quite small compared to 
the NLO case (respect to the leading order (LO) result), while another one is their dependence 
on the relevant scales (factorization, renormalization) of the process, which are arbitrary. The 
sensitivity on these scales is reduced by increasing the order of the perturbative expansion. 
A second point concerns the specific dependence of the result on the renormalization group 
evolution (RGE) and on the way we select the solution of the corresponding equations. When 
a generic equation is defined in terms of a power expansion in some parameter (coupling), its 
solution can be either power expanded in the same parameter or can be computed exactly, 
both approaches being legitimate options. This indetermination should be interpreted as a 
theoretical error which needs to be quantified. In fact, solving the RGEs in one way or in 
another is equivalent to a summation or to a resummation of the perturbative solution, with 
differences between the various methods which start appearing beyond a certain order. 

2 Brute-force and expansion-based solutions 

It is probably convenient, in order to understand the motivations for writing CANDIA, to briefly 
go over a classification of the ways the DGLAP equation is solved numerically and characterize 
the difference between "brute-force" and "expansion-based" solutions of these equations. 

A brute-force solution is obtained by simply discretizing the equation by a finite difference 
method, and neglecting all the issues of perturbative accuracy that we have just mentioned. 
Candia, on the other hand, is based on the implementation of analytical ansatze which have 
been shown [4] to be solutions of the evolution equations. These allow to keep track system- 
atically of the logarithmic corrections which are included in the final solution and to control 
directly its logarithmic accuracy. Being the logarithms accompanied by powers of a s , all the 
expanded solutions correspond to solutions of a given (and different) accuracy in the strong 
coupling. In particular, the program implements also resummed expansions, such as (1461) and 
(1601) . which involve logarithms of more complicated functions. 

Notice that in the DGLAP case - but the same issue appears in any equation whose right- 
hand-side is of a given accuracy in a parametric expansion - the kernel is only known up to 
a fixed order in a s , and therefore it is legitimate - in establishing the recursive form of the 
ansatz - to decide whether or not to drop all the information coming from the higher orders, 
orders over which we do not have complete control from the perturbative side and that are also 
present in the recursion relations. 

In this sense, it is clear that the "brute force" solution is just one of the many solutions 
which can be obtained by the analytical expansions. It can be reproduced from an expansion- 
based approach by extending the ansatz so to include all the logarithmic powers of the form 
a™ ln(a s /ao) n , with n, m integers. In practice n and m are finite integers, but chosen sufficiently 
large in the actual implementation of the numerical code. 

The approach followed in our program, though entirely formulated in x-space, shares some 
of the features which are typical of those Mellin methods that also rely on an analytical ansatz. 
In fact also in this case the exact solution is obtained recursively [5]. We also remark that 
the Mellin method, such as the one implemented in [1], and our method overlap as for overall 
treatment of the logarithms, though they implement different partial summations. The proof 
that the various analytical ansatze used by us are solutions of the evolution equations is obtained 
by construction, solving the recursion relations for the unknown coefficient functions of the 
expansion in Mellin space [6]. 

We have implemented two classes of solutions: the "exact" solutions and the "truncated" 
ones. The latter are obtained by retaining contributions to the expansion up to a certain power 
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of a s , which can be chosen by the user. We have found that the numerical change induced 
by the variation of the ansatz on some hadronic observables at NNLO is of the order of a 
percent. Though this variation is indeed small if compared to the change from LO to NLO of 
the cross section, it is comparable to its change from NLO to NNLO (quantified to be around 
3% on the Z peak [6]). For this reason, we think that Candia can be of help in the studies of 
resonant processes where a large amount of experimental data can be collected, such as on the 
Z peak. Here the small theoretical errors that are due either to the various ways of handling 
the evolution or to the changes induced by going from NLO to NNLO in the hard scatterings 
are far larger than the experimental statistical errors coming from the direct measurements at 
the LHC. In this work we are going to describe the basic features of our program and focus 
essentially on the evolution part. More tools which may be useful for QCD partonometry and 
for the search of extra Z' will be released in the near future. 



3 Summary of definitions and conventions 

We start by summarizing the notations and definitions that we will be using in the description 
of the program. The general mathematical structure of the DGLAP equation is 

f(x,v 2 F ) = P(x,a s (Li 2 F ))®f(x,Li 2 F ) (1) 



d In n F 

where the convolution product is defined by 



[a ® b] (x) = f % (*) b{y) = f ^a{y)b (*) (2) 

Jx y \yj Jx y \vj 

and the perturbative expansion of the kernels and of the beta function up to NNLO are respec- 
tively 

P(x,a s ) = (g) PW{ X ) + V>(a;) + G|) +0(at) (3) 



and 



where 



da s (n R ) _Po 2 _ Ji_ 3 _ Jh_ 4 5 
P[(Xs) dln/4 4tt s lGvr 2 ^ 64tt 3 s + [ s) [ ) 

A) = Y Nc ~t Tf (5) 

34 20 

Bi = -Nl - -N C T, - iC F T, (6) 
ft = 2 ^Ni + 2CiT,- 2 ^C F N c T l -l^NlT l + ^C F r i+ l -§N c T} (7) 

and 

N%-1 A 1 
Nc = S, C F = ^- = -, T f = T R n f = -n f . (8) 

Nc is the number of colors and n/ is the number of active flavors, selected by the mass condition 
m q < I^f, for a given factorization scale [i F . We also denote with \i R the renormalization scale. 
The evolution (DGLAP) kernels are distributions whose general form is given by 

P(x) = P^x) + + P 3 5(l - x), (9) 

U _ x )+ 
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with a regular part Pi(x), a "plus-distribution" part P2(x)/(1 — x) + and a delta-function term 
P 3 5(l — x). Given a continuous and differentiable function a(x) defined in the [0, 1) interval, 
but singular at x = 1, the action of the plus-distribution [a(x)] + is defined by 



f(x)[a(x)]+dx = / (f{x) -/(!)) a(x)dx, 



(10) 



where f(x) is a regular test function. Alternatively, an operative definition (that assumes full 
mathematical meaning only when integrated) is the following 



[a(x)] + = a(x) — 5(1 — x) / a(y)dy. 



(11) 



The Parton Distribution Functions (PDFs) appear in the evolution, in general, multiplied 
by a factor x. Eq. (DO) for f(x) = xf(x) then reads 







f(x, /4) = x [P(x, a s (n 2 F )) ® f(x, H 2 F )] . 



(12) 



We are now going to compute the convolution products that we actually evaluate in the program, 
having in mind the general form of the kernel Q. 
The treatment of the delta-function part is trivial 



x[P 3 5(l-x)®f(x)} =xP 3 
while for the regular part we get 



S(l-y)f 



Psf(x) 



(13) 



f^pA-)m=f^-p. (-)m. d4) 

Jx y \yj Jx y y \v) 



x [P\(x) ® f(x)\ = X 



Since the functions that we have to integrate are strongly varying or even singular at low x, to 
enhance the numerical accuracy we introduce a new integration variable z defined by y = x z , 
and Eq. ( TT4l) is mapped into 



x[P 1 (x)(S>f(x)] = -hyx / dzx 1 - z P 1 (x 1 ~ z )f(x 



(15) 



Finally, for the plus-distribution part we get, after some algebraic manipulations the relation 

P2(X 



X 



[1-x 

that with the mapping y = x z becomes 
P 2 (x 



1 (|// /M//)./V/ f ) /Mi )/(,■) + Ml)m ln(1 _ x) (16) 



X 



(1-x). 



■nx / ' d ^^>'"> - + ft(l)/(,) ln(l - x). (17) 







1 - x z 



The program, being entirely developed in (Bjorken's) x-space uses the iteration of the convo- 
lution in the form defined above both in the non-singlet and singlet sectors. 
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4 Non-singlet and singlet structure 



We start by defining the combinations 

qf ) = q i ±q i (18) 
n f 

^ = Eft (±) - ( 19 ) 

i=l 

From a mathematical point of view, the distributions belonging to the non-singlet sector evolve 
with a decoupled DGLAP equation of form (pQ), while the singlet combinations mix with the 
gluon distribution. The singlet DGLAP matrix equation is 



g ( <? (+) (^/4) \ = f P qq {x,a s {fi 2 F )) P qg (x,a s (fi 2 F )) \ / q {+ \x,fi 2 F ) 

din 11%, \ g{x,fi 2 F ) J \ Pgq{x, a s {fi 2 F )) P gg (x,a s (ii 2 F )) J \ g{x,/i F ), 

or, in vectorial notation 

U .2 \ Tt/ „ /,.2\\ „/„ ,.2 



(20) 



s(x, n F ) = P(x, a 8 (fi F )) ® s(x,n F ). (21) 



dm /if/ 

The general structure of the non-singlet splitting functions is given by 

P<li<lk = pQiQk = ^ikPqq + Pqq, (22) 
P<im = ^QiQk = ^ikPqq + Pqq- (23) 

where V and S are usually referred to as the "valence" and "sea" contributions. This leads to 
three types of non-singlet distributions which evolve independently: the flavor asymmetries 

(24) 

governed by the combinations 

PjVS = Pqq ^ Pqq ' (25) 

and the sum of the valence distributions of all flavors which evolves with 

PV S = P v qq - Pi + n f (P q s q - F£) = P NS + P* s . (26) 
Notice that the quark-quark splitting function P qq can be expressed as 

Pqq = PNS + n f (Pq S q + Pq'q) = P£ S + Pp.- (27) 

with ps denoting the so-called "pure singlet" terms. We remark that the non-singlet contribution 
is the most relevant one in Eq. ( 1271) at large x, where the pure singlet term P ps = P qq + P qq is 
very small. At small x, on the other hand, the latter contribution takes over, as xP ps does not 
vanish for x — > 0, unlike xP^ s . The gluon-quark and quark-gluon entries are given by 

Pq 9 = n f P m , (28) 

Pgq = P 9 q, (29) 
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in terms of the flavor-independent splitting functions P qi9 = P qi9 and P m = P ggi . With the 
exception of the first order part of P qg , neither of the quantities xP qg , xP gq and xP gg vanish for 
.r -> 0. 

In the expansion in powers of a s of the evolution equations, the flavor-diagonal (valence) 
quantity P q v q is of order a s , while P qq and the flavor-independent (sea) contributions P qq and 
P qq are of order a 2 s . A non- vanishing difference P qq — P qq is present at order af. 

The next step is to choose a proper basis of non-singlet distributions that allows us to 
reconstruct, through linear combinations, the distribution of each parton. The singlet evolution 
gives us 2 distributions, g and q^ + \ so we need to evolve 2n/ — 1 independent non-singlet 
distributions. At NNLO we choose 

1. q(~\ evolving with P^ s ; 

2 - QnIu = ~ ( for 2 < i < n f), evolving with P NS ; 

3 - Qns,u = <li +) ~ ( for 2 < z < n/), evolving with P^ s , 
and use simple relations such as 



to perform the reconstructions of the various flavors. Choosing i = 1 in (i30l ). we compute 
q[ ^ from the evolved non-singlets of type 1 and 2 and from the evolved singlet q^ and 
non-singlet of type 3. Then from the non-singlets 2 and 3 we compute respectively q\ ^ and 
g t - for each i such that 2 < i < rif, and finally qi and g^. 



and the same does each linear combination thereof, in particular q\ for each flavor i. The 
basis of the 2nf — 1 non-singlet distributions that we choose to evolve at NLO is 

1. q\~^ (for each % <«-/), evolving with -P^ 1 ^, 

2. g^5!j = g^ — gf 4 ^ (for each i such that 2 <i < rif), evolving with Ppcj 1 ^, 

and the same we do at LO, where we have in addition Pp^ = P N $ - , being P^q = 0. 

5 The algorithm 

We briefly review in this section the algorithm on which Candia is based. More details and a 
theoretical discussion can be found in our previous papers [4,6]. 

5.1 Non-singlet: exact solution 

The proof of equivalence between the logarithmic expansions implemented in Candia and the 
alternative approach based on the use of Mellin moments, as used in reference [1], is easily 




(30) 
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established at leading order. For this we recall the definition of the Mellin transform of a given 
function a(x), 

a(N) = [ a(x)x N ~ 1 dx, (31) 
Jo 

that maps convolution products into ordinary products 

[a®b](N) = a(N)b(N). (32) 

Traditional algorithms based on Mellin space solve the equations algebraically in N and then 
perform a numerical inversion using a saddle path evaluation of the complex contour. This 
technique is completely bypassed in our approach. One of the advantage of Candia is that a 
given initial condition for the PDFs, usually formulated in x-space, does not need to be fitted to 
a given functional form in moment space, which is instead typical of a given numerical imple- 
mentation of the Mellin algorithm. In fact, the functional form in moment space in some cases 
may even not be general enough, and may not allow the evolution of quite singular distributions 
at small x. From our experience, we have found that fitting special initial conditions in x-space 
forces the user, in codes based in Mellin space, to modify the evolution code by himself, with 
dubious results. Candia, by eliminating this unappealing feature of algorithms based in Mellin 
space, allows any initial condition to be considered and removes the initial numerical error due 
to the fit of the initial condition to the pre-assigned functional form in moment space. 

Having clarified this point, we introduce a single evolution scale Q, leaving to the next 
sections the discussion of the separation between the factorization and renormalization scales. 
Switching to a s as the independent variable, the DGLAP equation (Dp) is rewritten in the form 

df(N,a 8 ) P(N, 

a = —, Ti — \-f( N > a s)- ( 33 ) 

5.1.1 Leading order 

Inserting in Eq. (133]) the perturbative expansions of P(N,a s ) and (3{a s ) (Equations (jHJ and 
fll}) arrested at LO, we get 

" M - ^ P( ° ){N) M , ,34, 



da s " f° a 2 



which is solved by 



a. 



2p(0)(jy) 



f(N, a.) = f(N, a ) ( - ) * = f(N, a„) exp f_ 2P °^ ln S!t\ ( 35 ) 

where we have set ot s = a s (Q 2 ) and «o = ®s(Qo)- Performing a Taylor expansion of the 
exponential we get 

The x-space logarithmic ansatz, that parallels this solution is 

t n (x) a s (Q 2 
n\ oc s {Ql) 



n=0 
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where the A n (x)'s are unknown functions. Setting Q = Q in (1371) we get the initial condition 

f(x,Ql)=A (x). (38) 

Inserting our ansatz (1371) into the DGLAP equation (0Q) and using the expansion of the kernels 
and of the beta function (J3J [4J arrested at the first term, after some algebra we derive the 
recursion relation 

A n+1 (x) = -^-[P^®A n ](x). (39) 
Po 

In the code, the value of A n (xk) = xA n (xt) for the PDF with index i (see TableE]), is stored 
in the variable A [i] [n] [k] ; x^ is the x-grid point stored in X [k] . 

5.1.2 Next-to-leading order 

Let's now move to the higher orders. 
At NLO Eq. {M} reads 

df(N,a s ) _ fe)PW(iV) + fe) 2 P (1) (AQ 

s 4f s ' I6vr 2 s 



the solution of which is 

'a s \ f Anp + a s pA ft 



f(N,a s ) = f(N,a )( 



2 p(0)(jv) 2p(°>(JV) 4pW(N) 



where we have set 



\a J \4tt/3o + a Q /3i 

/(A^o)e aW V^ M 

\n=0 / \m=0 



L = In— (42) 
«o 

M = ln f^ + a f (43) 

a(AT) = J^J (44) 

c(iV) = (45) 

Po Pi 



We then assume an x-space solution of the form 



,ra=0 ' / \m=0 



sr^ sr^ A n (x) C s - n (x) , 
2-^2-^ n\ (s-n)\ 

s=0 n=0 v ' 

s=0 n=0 v ' 



where in the first step we have re-arranged the product of the two series into a single series 
with a total exponent s = n + m, and in the last step we have introduced the functions 



B s n (x) = A n (x)C s . n (x), (n < s). (47) 

Setting Q = Q in ( 1461 ) we derive the initial condition on the recursive coefficients 

f(x,Q 2 ) = B° (x). (48) 

Inserting the ansatz (l4*6l) into the DGLAP equation ([TJ), using the expansions (J3], SJ) arrested at 
the second terms, and equating the coefficients of a and a 2 , we find the recursion relations 

B s n (x) = ~ [P<°> ® B s n z\] (x) (49) 
Po 

B° n (x) = -B s n+l {x)-^[P^®B s - l ]{x). (50) 

Pi 

These relations allow to compute all the coefficients B^ (n < s) up to a chosen s starting from 
Bq, the value of which is given by the initial conditions. Eq. (l4*9l allows to follow a diagonal 
arrow in the scheme shown in Table [U Eq. (l50l instead allows to compute a coefficient once 
we know the coefficients at its right and over it (horizontal and vertical arrows). If more than 
one recursion relation can be used to compute a coefficient B^, the program follows the fastest 
recursion chain to determine these, i.e. (I49l that involves instead of ?W. For each s the 
code does the following: 



1. computes all the coefficients B s n with n^O using ((19 

2. computes the coefficient Bq using (loOl) . 



In the code, the value of B^(xk) = xB^(xk) for the PDF with index i (see Table [5|), is 
stored in the variable B [i] [s] [n] [k] . 

B° 

I \ 

B l *~ B \ 

I \ \ 

r2 r2 r2 

B q <— B x k> 2 



jd3 , jd3 jd3 jd3 

B Q B 1 B 2 B 



Table 1: Schematic representation of the procedure followed to compute each coefficient B, 



5.1.3 Next-to-next-to-leading order 

At NNLO Eq. ((331) reads 

df(N,a s ) (^)P (0) (N) + (^) 2 PW(N) + {^) 3 P^(N) 
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f(N,a s ), (51) 



the solution of which is 



f(N, a s ) = f(N,a )e a ^ L e b ^ M e d ^ T 



\n=0 ' / \m=0 ' / \p=0 F ' / 

where we have introduced the definitions 

L = In— (53) 

_ 167T 2 /3 + + a% 

16tt 2 /3 + 4vra p\ + og/? 2 

^ 1 , 27rK-a )VW2-/3i 

— — = arctan — — — ; — — — — 55 

^A(3q(3 2 - (51 2n(8n(3 + (a s + a )p\) + a s a /3 2 

«(N) = -« ,50, 
Po 

x P(°)(JV) 4P( 2 )(iV) / 

6(iV) = _ *U. (57) 

Po P2 

d(JV) = ^lp(°)(iV)-8P«(iV) + ^P (2) (iV). (58) 

PO P2 

Notice that, if 4/3 /52 — /?i < (it occurs for rif = 6), T has to be analytically continued 



T = 1 octant 

y/Pf - 4/? /3 2 27r(87T/3 + (a 8 + a )A) + a s a /3 2 

We then assume an x-space solution of the form 

\n=0 ' / \m=0 ' / \p=0 ^' 

y^y^ A- t (x) gt-n(g) g»-t(g) L n M t-n T s-t 

2-^2^2^ n \ (t- n )\(s-t)\ 

s=0 t=0 n=0 v ' v 7 

oo s £ / \ 

= YYY u w .. L n M t ~~ n T s ~ t , (60) 
n\(t -n)\(s -t)\ ' 

where in the first step we have transformed the product of three series into a single series in the 
total exponent s = n + m + p, and we have set t = n + m. In the last step we have introduced 
the functions 

CUx) = A n (x)B t „ n (x)D s _ t (x), (n<t< s). (61) 



Setting Q = Q in ( 1601) we get the initial condition 

f{x,Ql)=C%{x). (62) 

Inserting the ansatz ([601 into the DGLAP equation (0Q) and using the expansions of the kernel 
and the beta function (El HI) arrested at the third order, equating the coefficients of a, a 2 and 
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a 3 we find the recursion relations 

CUx) = ~ [P<°) ® CJT^J (x) (63) 
Po 

= l C UW-^[ P(2) ^Hn]W (64) 

q n (x) = -2/3 1 (^ +1) Jx) + ^ +1)n+1 (a;))-8[p( 1 )®q; 1 ](x). (65) 

Also in this before, when we have to compute a given coefficient Cf n , if more than one 

recursion relation is applicable, the program follows the fastest recursion chain, i.e. in order 
(l63l . (16*51) and ([""""") ■ At fixed s the algorithm performs the following steps: 

1. computes all the coefficients C t s n with n ^ using ( l63l) : 

2. computes the coefficient C| using ("64") : 

3. computes the coefficient C t s with f ^ s using ( T65i ). in decreasing order of t. 

This procedure is exemplified in the scheme shown in Table [2] for s = 4. 

In the numerical program, the value of C^ n (xi~) = xCf n (xk) for the PDF with index i (see 
Table [""J), is stored in the variable C[i] [s] [t] [n] [k] . 









°4,4 






°3,3 


/^y4 
U A,3 




°2,2 


/^4 

°3,2 


r~iA 
°4,2 


1/ 


/~iA 
b 2,l 


/nr4 
°3,1 


/nr4 
' i 


riA , 
°1,0 ^~ 


~~ °2,0 


~~ °3,0 


riA 
~ °4,0 



Table 2: Schematic representation of the procedure followed to compute each coefficient C\ 
for s = 4. The underlined coefficients are computed via Eq. ( 1631) . 



5.2 The truncated solutions 

Besides the class of solutions of Eq. (pp) described in Section 15.11 which we have called exact, 
there is another important class of solutions, that we will call truncated, which are interesting 
since they correspond to x-space solutions which are accurate up to a certain order in a s . 
While in an exact solution all the logarithmic structures are resummed into few functions 
(Equations (142143ft at NLO and (I53II54II55I) at NNLO), the truncated ones are characterized by 
expansions in terms of simple logarithms of a s /ao, retained up to a chosen order. We give below 
some details on these expansions. Notice that all the recursive solutions built in the singlet 
sector - except for the lowest order ones, which are exact in any approach, either in Mellin 
space or in x-space - are of this type. In this sector we can build solutions of increased accuracy 
by using truncated solutions of higher order. To briefly discuss these types of solutions, we 



12 



perform an expansion in ot s of the quantity Pj (3 in the Eq. (1331) and re-arrange the DGLAP 
equation (NLO or NNLO) into the form 



di{N,a a ) 1 



— [Rq + a s Ri + a 2 s R 2 + ■■■ + a£Rj f(N, a,), 

(66) 



where we have considered the expansion up to order k and we have defined the following linear 
combinations of the P(iV) kernels 

Ro = - J-P (0) 

Po 

R --2^^ P(1) - P(0) ^] 

1 /PW Rip\ RoP 2 \ 
R2 = -- ( ^ + ^ + T6^J • 

(67) 



Ti \2ttPo 4/3 16tt/? c 



We call Eq. (l66l) the truncated version of the DGLAP equation, both in the singlet and non- 
singlet cases. 

The truncated equation, in Mellin space, can be solved in closed form, at least in the non- 
singlet case, obtaining a solution which is different from the exact solution of Eq. ( !33l) discussed 
before and having the following general form 

f(N, a,) = (^j exp |Ri(« s - a ) + ^R 2 («' -«') + ... + ^R K « - a£) | f(N, a ) . 

(68) 

Even after a truncation of the equation, the corresponding solution is still affected by higher 
order terms, and can be truncated. To obtain the truncated version of this solution is then 
necessary a further Taylor expansion around the point (a s , ato) = (0,0) for the two couplings 
-initial and final- 

R ° r 11 i 

f{N,a ). (69) 



1 + Ri(a s - a ) + ^R 2 (c^ - «o) + • • • + - a o, 

Z K 



Again, this expression holds in both singlet and non-singlet cases, thus we can generate truncated 
solutions for the parton densities. 

Passing to the singlet case, which is more involved, the truncated solutions of Eq. (1661 ) in 
Mellin space can be obtained (for a review see [1,4,6]) by the use of the [/-matrix ansatz. 
Basically this method consists of an expansion of the general solution around the LO solution 
as 

f(N, a s ) = [1 + a s V 1 (N) + a 2 s U 2 (N) + ... + a K 3 V K (N)] L{a s , a , N) 
[1 + a U 1 (N) + a 2 U 2 (N) + ... + « K U K (iV)] _1 f (JV, a„), 

(70) 

where the L(a s ,ao,A^) is defined as 

L(a s ,a ,N)= (^-) (71) 
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with the U operators defined through a chain of recursion relations, obtained by substituting 
Eq. (1701) into the truncated equation 



[R ,V 1 ] =U 1 -R 1 , 

[R ,U 2 ] = -R 2 -R 1 -U 1 + 2U 2 

; (72) 



With a Taylor expansion of the second line in Eq. (|70|) we obtain the truncated solution of order 
k for the singlet/non-singlet case in Mellin space. Once we have obtained such a solution, the 
x-space result is achieved by a Mellin inversion. This and the x-space approach merge as we 
increase the order of the expansion. Working directly in x-space, one can generate truncated 
solutions up to certain order k in a very general way by considering the following logarithmic 
ansatz 



i=0 



... (?3) 

n=0 k L*=0 J K ^ 

which is stunningly simple. One can demonstrate [4, 6] that inserting this expression in the x- 
space version of Eq. ( 1661 ) the solution obtained is equivalent to that of the £7-matrix prescription, 
in the sense that the terms of the two expansions are the same as far as the two expansions are 
implemented up to a sufficiently large order. In our approach the solution is expanded in terms 
of logarithms of a s /a and one controls the accuracy by inserting powers of ot s . The coefficients 
S l n (x) are determined by solving the chain of recursion relations generated by inserting the 
logarithmic ansatz in the DGLAP, grouping the terms proportional to the same power of a s 
and neglecting the terms in and higher powers. 

We obtain the following recursion relation for the S° coefficient 

S° +1 (*) = -|-[P (0) ®S°](x), (74) 

Po 



while for S 1 we have 



~ [P (0) ® Si] (x) - \ [P« ® S°] (x) . (75) 
The recursion relations for the S* coefficients (i = 2, 3, . . . , ac) are, at NLO 



P 1 ai-l /„\ .-o» r„\ a i \ @- 



^n+i( x ) — ~ A 7rSn+i( x ) iS l n (x) (i 1)- -r-SJ, 1 (x 
4vrp 47rPo 



A [P<°> ® Si] (x) - -L [PW ® S- 1 ] (x) (76) 



and at NNLO 



s " +i(x) 4^ s "+ i(a;) i^ s "+ i(a:) 



-isi(x) - (,• - ij-^-sjr^x) - (i - 2 )-^-s*r a (x) 

-A [P<°> 6) Sy (x) - 4" [P (1) ® Sr 1 ] (*) - ^TT [ p(2) ® S tT 2 ] (*) (77) 

PO T^PO - Z7T Po 
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These relations hold both in the non-singlet and singlet cases and can be solved either in x- 
space or in iV-space in terms of the initial conditions f (x, «o) — A further check of the 

overlap of the two methods is obtained numerically. 

Since there is no closed form solution for the singlet DGLAP equation, we always use the 
method described in the current section in this specific sector. On the other hand, we leave to 
the user the possibility to choose the method of solution for the non-singlet equations, exact 
or truncated that they can be, and this choice can be made at compilation time (see Section 
I8.2.2H . If one chooses the exact method implemented in candia.c, we need to define the S 
coefficients of eq. (173]) only for the singlet sector: 

S n ,singlet = ( Qi ^ j • (78) 

In the code, the value of S^Xk) = xS l n (xk) for the gluon is stored in the variable S [i] [0] [n] [k] 
and for in S[i] [1] [n] [k] . If one chooses instead to solve both the singlet and the non- 
singlet with the truncated method (implemented in candia_trunc . c) the array S has to ac- 
commodate also the non-singlet distributions, and the value of S % n (xk) for the PDF with index 
j (see Table E|) is stored in the variable S[i] [j] [n] [k] . 

We remark that in a previous NLO version of the implementation of the algorithm both the 
singlet and the non-singlet sectors had been solved using truncated ansatz with k = 1 [7]. 



6 Renormalization scale dependence 



As we move to higher orders in the expansion of the kernels, the presence of the strong cou- 
pling constant a s allows an independent renormalization scale (jlr on the right-hand-side of the 
evolution equation. This dependence is, in general, completely unrelated to the factorization 
scale. Thus, we can formally rewrite the DGLAP equation as follows 







<91n ftp 



f(x, fi%, /4) = p{x, /4, /4) ® /(^ /4> /4) 



(79) 



where the splitting functions have acquired a /i# dependence simply by expanding a s (fi 2 F ) in 
terms of a s (fJ, R ) 



a s (fx 2 R ) = a s (/40 



-<wf)-i— + Ta \9 \~Po L - Pl L ) 



Air 



(47T) 



(80) 



The explicit replacements of the kernels in this new re-organization of the defining equation are 
given by 



p(0)( x ) ^ p(0)( x ) 

pW( x ) _> p<X)( x ) 



^P^{x)L 



p(2)( x ) _p(2)( x ) -P Q LPW{X) 

where the logarithmic structures are identified by 







(81) 



L = ln4- 



(82) 
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Candia allows a determination of the evolution with the two scales held separate //p ^ I^r- In 
particular we have implemented the case in which fip and are proportional. The proportion- 
ality factor kn = ^%/^ R is entered by the user as a command-line argument (see Section [8.2.3j) . 
One can easily figure out, practicing with the program, that on increasing the perturbative 
order the dependence on the renormalization/factorization scale reduces. 



7 Treatment of the quark mass thresholds 

If Candia is executed with the variable flavor number scheme (VFNS) option (fns set to 1 
in the command line) and with the macro HFT set to 1 (that is the default value) in the file 
constants. h, the program will implement the matching conditions with rif and rif + 1 light 
flavors both for the running coupling a s and for the parton distributions. The transition from 
the effective theory with rif light flavors to the one with rif + 1 is made when the factorization 
scale reaches the renormalized pole mass of a heavy quark, fj,p — m nf+ i. The matching condition 
up to NNLO for the running coupling is [8] 



or, otherwise 

(n / +l) _ (rif) 



2 L 
6n 



3 f 1 



7T 



as " + 



a 



(n f ) 



I 2 L 
6tt 



(«/)l 3 f 1 



48tt 2 



— - —L — 

36 ~ 24 ~ 24 



14 + 38L 



O 



a 



("/) 



O 



a 



(n f ) 



1 1 



where L = \nk R = ln(//|j/ n 2 F ) . 

The matching conditions for the parton distributions up to NNLO are [9] 



(83) 



(84) 



a 



(»V+i) ' 



47T 



'A N qq S f®lt f) 



(x) 



(85) 



where l=q, q and i = 1,2, ... ,Uf, are the light quark/antiquark flavours; 



a 



(n/+l) ' 



4:71 



A*™ (S 



(x) + 



(86) 



for the gluons, and 



$' +1) (*) = g£7 +1) (s) 



(n/+l)' 



47T 



^(2) ^ q(+Unf) 



(x) + 



IJf®^] Or)} (87) 



for the heavy flavors. 



8 Description of the program 

8.1 Content of the package 

The code is unpacked with the command 
tar zxvf candia_l . . tar .gz 
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that will create the directory candia_1.0, containing the following files. 

• candia. c and candia_trunc . c are the files including the main function, each one imple- 
menting a different method of solution: the former solves the non-singlet sector with the 
exact solution, while the latter uses the truncated method. 

• makefile and makef ile_trunc are the corresponding makefiles. If one is using other 
compilers than gcc and gf ortran these files need to be edited. 

• constants . h is a header file containing some parameters that the user may want to edit. 

• xpns2p.f and xpij2p.f are Fortran codes [10,11] in which a parametrized form of the 
NNLO kernels is defined. Very few modifications have been done to make them compatible 
with our code. 

• hplog.f is a Fortran code [12] in which a subroutine that computes numerically the 
harmonic polylogarithms up to weight 4 is implemented. Harmonic polylogarithms are 
defined in [13]. 

• partonww.f is just a merging of the three Fortran codes mrst20011o .f , mrst2001.f and 
mrstnnlo.f by the MRST group [14,15] to access their grids of LO, NLO and NNLO 
parton densities. Very few modifications have been done. 

• lo2002.dat, alfll9.dat and vnvalfll55.dat are the MRST parton densities grids at 
LO, NLO and NNLO respectively. 

• a02m.f is the Fortan code by Alekhin [16] to access his grids of LO, NLO and NNLO 
parton densities. 

• a02m.pdf s_i_vfn and a02m.dpdf s_i_vfn with % — 1, 2, 3 are the files in which the grids 
of the Alekhin parton densities in the variable flavor number scheme are stored. 

8.2 How to run the program 

Let us describe the different steps that the user will encounter in a run of Candia. 



8.2.1 Editing the file constants. h 

The header file constants. h contains some macros and two arrays that the user may want to 
change before compiling. The macros are described in Table [3] 



Macro 


Default 


Description 


GRID_PTS 


501 


Number of points in the x-grid 


NGP 


30 


Number of Gaussian points 


ITERATIONS 


15 


Number of iterations 


INTERP.PTS 


4 


Parameter used in the polvnomial interpolation (see Section 18. 3. 3D 


HFT 


1 


Switch for the heavy flavors treatment [8,9] 



Table 3: Macros defined in constants. h. 
Besides these macros, two important arrays are defined in constants. h. 
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• xtab. The values of x listed in this array need to be in increasing order. The first (lower) 
value is the lower value of x in the grid. The last (upper) value must be 1. The eventual 
intermediate values will be forced to be in the grid. This array must contain at least two 
values. 

• Qtab. An output file with the PDF values computed at the end of the evolution will be 
generated for all the values of listed in this array. 

8.2.2 Compiling 

At this point the user has to choose the solution method for the non-singlet sector. 

The main file in which the exact method is implemented and the corresponding makefile 
are candia.c and makefile, while the truncated method is implemented in candia_trunc . c 
and the corresponding makefile is makef ile_trunc. 

In Table [4] we show for each method of solution the compilation command that the user 
should type and the executable file that will be produced. 



Method 


Command 


Executable 


exact 


make 


candia.x 


truncated 


make -f makef ile_trunc 


candia_trunc . x 



Table 4: Summary of compilation commands. 



8.2.3 Running 

Let us suppose that the user wants to use the exact method of solution, so the executable 
is called candia.x. The whole procedure applies for the truncated method as well, replacing 
candia.x with candia_trunc . x. 

The user has to supply six parameters to the command line. To have a quick usage update 
he/she can just type the line 

candia.x 

whose output is self-explanatory. 
USAGE 

candia.x <perturbative_order> <truncation_index> <input_model> <kr> <fns> <ext> 
<perturbative_order> can be 0, 1 or 2 

<truncation_index> cannot be less than <perturbative_order> 
<input_model> : 

0= Les Houches toy model 

1= MRST parametrization 

2= MRST grid at 1.25 GeV~2 (minimal value in their grid) 

3= Alekhin parametrization 

4= Alekhin grid at 1 GeV~2 
<kr> is the ratio mu_r~2 / mu_f~2 
<fns>: 
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0= fixed flavor number scheme 
1= variable flavor number scheme 
<ext> is the extension of the output files (max 3 characters allowed) 

For example, to run the program with the exact solution method at NNLO, with a truncation 
index k — 6 and the MRST grid as input, with /i R — /i 2 F , in the variable flavor number scheme 
and choosing the extension dat for the output files, one should type the command 

candia.x 2 6 2 11 dat 

8.2.4 Understanding the output files 

For example, if we choose dat as an extension of the output files and our Qtab array is defined 
by the following line in constants. h 

const double Qtab [] ={50 . , 100 . , 150 . , 200 . } 

when the program exits, the user will find in the directory some files called Qi.dat, with i 
ranging from to the number of elements of Qtab (4 in our example). Q0.dat is the summary 
of the relevant output files, and in our example it will look like 



file Q Q~2 



Qi.dat 50 2500 

Q2.dat 100 10000 

Q3.dat 150 22500 

Q4.dat 200 40000 



The other files (for example Qi.dat, showing the PDFs at 50 GeV) have one line for each 
point in the x-grid and 14 columns: x, xg, xu, xd, xs, xc, xb, xt, xu, xd, xs, xc, xb and xi. 
To print the output in a different form, one should edit a small portion of the code. 

8.3 Functions 

8.3.1 Safe allocation functions 

void *Malloc (size_t size); 

void *Calloc (size_t nmemb,size_t size); 

These functions have the same syntax as malloc and calloc of Standard C. The only 
difference is that they check if there is enough memory available to perform the allocation, and 
if not they terminate the execution of the program. 

8.3.2 Function gauleg 

void gauleg(double xl, double x2, double *x, double *w,int n) ; 

This function is taken from [17] with just minor changes. Given the lower and upper limits 
of integration xl and x2, and given n, gauleg returns arrays x[0 , . . . ,n-l] and w[0 , . . . ,n-l] 
of length n, containing the abscissas and weights of the Gauss-Legendre n-point quadrature 
formula. 
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8.3.3 Interpolation functions 
double interp (double *A, double x) ; 

double polint (double *xa,double *ya,int n,double x) ; 

Given an array A, representing a function known only at the grid points, and a number x in 
the interval [0, 1], interp returns the polynomial interpolation of grade INTERP_PTS * 2 — 1 of 
the function at the point x through an appropriate call of polint. The number of points used 
for the interpolation are controlled by the parameter INTERP_PTS. This identifies the number of 
grid points preceding and following x which are globally used for the interpolation (for a total 
number of IMTERP_PTS * 2). If the number of available grid points before or after x is smaller 
than the value INTERP_PTS, more points at smaller or larger values of x are chosen, to ensure 
that the total number of used points is always INTERP_PTS *2. The actual interpolation is done 
by polint, a slightly modified version of the function with the same name that appears in [17]. 

8.3.4 Function convolution 

double convolution(int i,double kernel (int , double) , double *A) ; 

Given an integer i, to which corresponds a grid point x«, a two variable function kernel (i , x) , 
representing a kernel P(x), and an array A, representing a function xf(x) = f(x) known at 
the grid points, convolution returns x[P ® f](x) as the sum of the three pieces (fT5l) . (fT71) and 
(TT31) . The integrals are computed using the Gauss-Legendre technique. 

8.3.5 Function RecRel.A 

double RecRel_A (double *A,int k, double P0 (int , double) ) ; 

Given an array A, representing a function A(x) known at the grid points, an integer k (to 
which corresponds a grid point x k ) and a two variable function P0(i ,x) , representing a leading 
order kernel, RecRel_A returns 

- A [P(°> ® A] (x k ) (88) 
Po 

i.e. the RHS of Eq. (1391 ) (or equivalentlv the RHS of each component of the vectorial equation 

(3D)- 

8.3.6 Function RecRel.B 

double RecRel_B (double *A, double *B,int k, 

double PO (int , double) , double PI (int , double) ) ; 

Given two arrays A and B, representing the functions A(x) and B(x) known at the grid 
points, an integer k (to which corresponds a grid point x k ) and two functions P0(i,x) and 
Pl(i,x), representing respectively the LO and NLO part of a kernel, RecRel_B returns 

- 1- [?(") ® b] ( Xk ) - 4r [ p{1) ® A ] M (89) 

Po ttPo 

i.e. the structure of the convolution products that appear in Eqs. (1751) and ( 1761 ). 

When 7^ /j,f, some additional terms are included to take into account the transformation 
rules of the kernels (|8TI) . 
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8.3.7 Function RecRel.C 

double RecRel_C (double *A, double *B,double *C,int k,double PO (int , double ) , 
double PI (int , double) , double P2 (int , double) ) ; 

Given three arrays A, B and C representing the functions A(x), B(x) and C(x) known at the 
grid points, an integer k (to which corresponds a grid point Xk) and three functions P0(i,x), 
Pl(i,x) and P2(i,x), representing respectively the LO, NLO and NNLO part of a kernel, 
RecRel_C returns 

- J" [^ (0) ® C] (x k ) - -L [PW ® B] (x k ) - -L- [P^ ® A] (x k ) (90) 
Po ' ^Po ' ^ Po 

i.e. the structure of the convolution products that appear in Eq. (1771) . 

When ^ fip, some additional terms are included to take into account the transformation 
rules of the kernels (1ST1 ). 

8.3.8 Function RecRel_Diag 

double RecRel_Diag (double *C,int k,double PO (int , double) ) ; 
The same as RecRel_A. 

8.3.9 Function RecRel.Vertl 

double RecRel_Vertl (double *B,int k,double PO (int ,double) ,double PI (int ,double) ) ; 

Given an array B, representing a function B(x) known at the grid points, an integer k 
(to which corresponds a grid point x k ) and two functions P0(i,x) and Pl(i,x), representing 
respectively the LO and NLO part of a kernel, RecRel_Vertl returns 

- i- [P« 9 B] (x k ) (91) 
Pi 

i.e. the convolution product that appears in Eq. (i50l) . 

When /ir/ /J,p, some additional terms are included to take into account the transformation 
rules of the kernels (I5T1) . 

8.3.10 Function RecRel_Vert2 

double RecRel_Vert2 (double *C,int k,double PO (int , double ) , 

double PI (int , double) , double P2 (int , double) ) ; 

Given an array C representing a function C(x) known at the grid points, an integer k 
(to which corresponds a grid point x k ) and three functions P0(i,x), Pl(i,x) and P2(i,x), 
representing respectively the LO, NLO and NNLO part of a kernel, RecRel_Vert2 returns 

- ^ [P (2) S> C] (x k ) (92) 

P2 

i.e. the convolution product that appears in Eq. (i64l) . 

When /Ir ^ /j,f, some additional terms are included to take into account the transformation 
rules of the kernels (|B"T1) . 
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8.3.11 Function RecRel_Horiz 

double RecRel_Horiz (double *C,int k,double PO (int ,double) , double PI (int ,double) ) ; 



Given an array C, representing a function C(x) known at the grid points, an integer k 
(to which corresponds a grid point Xk) and two functions P0(i,x) and Pl(i,x), representing 
respectively the LO and NLO part of a kernel, RecRel_Horiz returns 

- 8 [P (1) ® C] (x k ) (93) 



i.e. the convolution product that appears in Eq. (I6al) . 

When Hr^ 1 Hf-> some additional terms are included to take into account the transformation 
rules of the kernels (I5T1) . 



8.3.12 function 

double BetaO(int f ) ; 
double BetaKint f ) ; 
double Beta2(int f ) ; 
double Beta(int order, double alpha); 



Given an integer f , representing the number of active flavors n/, BetaO, Betal and Beta2 
return respectively (3q, (3\ and @2 as in Eqs. (jHj), ((6)) and ([7]). 

Given the perturbative order and the value alpha of a s , Beta returns the value of the (3 
function given by Eq. (jlj) arrested at the chosen order. 

8.3.13 Function alpha_rk 

double alpha_rk(int order, double alphaO , double Qi, double Qf ) ; 



Given the perturbative order, an energy scale Qi (Qi) where a s has the known value alphaO 
(qjo) and another energy scale Qf (Qf), alpha_rk returns ce s (Qf), i.e. the solution of the Cauchy 
problem 

^ =«<«.). («) 

where t = lia(Q 2 ), with the boundary condition a s (t ) = a s (hx(Q 2 )) = olq. The differential 
equation has a closed form solution only at LO; at higher orders, a fourth order Runge-Kutta 
method is used. To get a correct result, the number of active flavors at both energy scales must 
be the same. 



8.3.14 Functions post and pre 

double post (double pre,int order); 
double pre (double post,int order); 

These functions implement the matching condition at a quark mass threshold for a s , as in 
Eqs. ( 1841 ) and ( 1831 ) arrested at the given perturbative order, post returns as in Eq. (184D 

given (pre), while pre returns ai n/ ^ as in Eq. (1831) given a^ 1 ^ (post). 
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8.3.15 Special functions 



double Li2 (double x) 
double Li3 (double x) 
double S12 (double x) 



These functions return respectively 



Li 2 {x) 



ln(l-y) 



dy 



(95) 



y 



Li 3 (x) 




(96) 




/ 

Jo 



\n 2 (l-y) 



dy 



(97) 



y 



using the program hplog [12]. 

8.3.16 Function fact 
double fact(int n) ; 

This function returns the factorial n\ 

8.3.17 MRST initial distributions 

double xuv_l(int order, double x) ; 
double xdv_l(int order, double x) ; 
double xS_l(int order, double x) ; 
double xg_l(int order, double x) ; 
double xD_l(int order, double x) ; 

Given the perturbative order and the Bjorken variable x, these functions return respectively 
xuy(x), xdv(x), xS(x), xg(x) and xA(x) in parametric form at Ql — 1 GeV 2 as defined in [15] 
and [14]. 

8.3.18 Alekhin initial distributions 

double xuv_2(int order, double x) ; 
double xus_2(int order, double x) ; 
double xdv_2(int order, double x) ; 
double xds_2(int order, double x) ; 
double xss_2(int order, double x) ; 
double xg_2(int order, double x) ; 

Given the perturbative order and the Bjorken variable x, these functions return respectively 
xuv(x), xu 3 (x), xdv(x), xds(x), xs s (x) and xG(x) in parametric form at Ql = 9 GeV 2 as 
defined in [16]. 
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8.3.19 Initial distributions for the Les Houches toy model 

double xuv_3 (double x) ; 
double xdv_3 (double x) ; 
double xg_3 (double x) ; 
double xdb_3 (double x) ; 
double xub_3 (double x) ; 
double xs_3 (double x) ; 

Given the Bjorken variable x, these functions return respectively ), xg(x), 

xd(x), xu(x) and xs(x) in parametric form at Ql = 2 GeV 2 as defined in [18]. 

8.3.20 LO kernels 

double P0MS(int i, double x) ; 

double P0qq(int i, double x) ; 

double P0qg(int i,double x) ; 

double P0gq(int i, double x) ; 

double P0gg(int i, double x) ; 

Given the Bjorken variable x, these functions return, depending on the value of the index i 

1. the regular part of the kernel Pi(x); 

2. the plus-distribution part of the kernel p2(x); 

3. the delta-function part of the kernel P 3 , 

as in Eq. ([9]). The kernels returned are respectively P^g, Pgq\ Pqg\ Pgq 

and P$ [19]. 

8.3.21 NLO kernels 

double PlNSm(int i,double x) ; 
double PlNSp(int i, double x) ; 
double Plqq(int i, double x) ; 
double Plqg(int i, double x) ; 
double Plgq(int i, double x) ; 
double Plgg(int i,double x) ; 

Given the Bjorken variable x and the index i with the same meaning of Section l8.3.20l these 
functions return respectively P^ , Ppg^, Pqq ', Pqg , Pgc? and Pgg [19]. 

8.3.22 NNLO kernels 

double P2NSm(int i, double x) ; 
double P2NSp(int i, double x) ; 
double P2NSv(int i, double x) ; 
double P2qq(int i, double x) ; 
double P2qg(int i,double x) ; 
double P2gq(int i, double x) ; 
double P2gg(int i, double x) ; 

Given the Bjorken variable x and the index i with the same meaning of Section l8.3.20l these 
functions return respectively P^s ■> Pps i ^NS j Pqq ■> Pqg \ Pgq an d Pgg m a parametrized 
form [10,11], using the Fortran files xpns2p.f and xpij2p.f. 
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8.3.23 Quark mass threshold kernels 

double A2ns(int i, double z) ; 

double A2gq(int i, double z) ; 

double A2gg(int i,double z) ; 

double A2hq(int i, double z) ; 

double A2hg(int i,double z) ; 

Given the Bjorken variable z and the index i with the same meaning of Section l8.3.20l these 
functions return respectively A^ 2 \ A^, A^ ( jJ, i§J 2) and A^f [9]. 

8.4 Distribution indices and identifiers 

In Table [5] we show the numerical index associated to each PDF appearing in the program. To 
the PDF labelled by the number index we associate a string id [index] , reported in the third 
column. Although not used in the current implementation, these identifiers can be useful if 
one wish to modify the output of the program: for example if one wants to print a file for each 
PDF, these strings can be used as part of the file name. 






gluons, g 


g 


1-6 


quarks, q i7 sorted by increasing mass (u, d, s, c, b, t) 


u,d,s,c,b,t 


7-12 


antiquarks, q { 


au , ad , as , ac , ab , at 


13-18 


it' 


urn , dm , sm , cm , bm , tm 


19-24 


#> 


up,dp,sp,cp,bp,tp 


25 


„(-> 


qm 


26-30 


dd, sd , cd,bd,td 


31 


q (+) 


qp 


32-36 




ds , ss , cs ,bs ,ts 



Table 5: Correspondence between indices, parton distributions and identifiers. 



8.5 Main variables 

In Table [6] we describe the main variables defined in the program. 

9 Results 

In Tables [7] and M we compare the results obtained running Candia with the exact and the 
truncated method, respectively, with those obtained using Pegasus [1] in the variable flavor 
number scheme and at a final scale \ip = fiR = 100 GeV. The initial conditions are taken 
from the Les Houches toy model with Ql = 2 GeV 2 [18]. We have made these runs using the 
following constants. h file 

#define GRID.PTS 801 
#define WGP 30 
#define ITERATIONS 15 
#define IMTERP_PTS 4 
#define HFT 1 
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X[i] i-th grid point, Xi 

XG[i] i-th Gaussian abscissa in the range [0, 1], Xi 

WG[i] i-th Gaussian weight in the range [0, 1], Wj 

order perturbative order 

trunc truncation index k 

input identifier of the input model (see Section I8.2.3H 

kr k R = h 2 r /h 2 f 

ntab[i] position in the array X of the x value xtab[i] (i.e. X [ntab [i] ] =xtab [i] ) 

nf , Nf number of active flavors, rij 

nf i number of active flavors at the input scale 

nf f number of active flavors in the last evolution step 

betaO Pq 

betal (3i 

beta2 (3 2 

beta (3 

log_mf2_mr2 In fi 2 F /fi R 

Q[i] values of fip between which an evolution step is performed 

alpha_pre [i] value of a s just below the i-th flavor threshold 

alpha_post [i] value of a s just above the i-th flavor threshold 

alphal a s (/i^), where \i % R is the lower Hr of the evolution step 

alpha2 a s (n R in ), where fi R n is the higher [i R of the evolution step 

A[i] [n] [k] coefficient A n {x^) for the distribution with index i 

B[i] [s] [n] [k] coefficient B s n {xk) for the distribution with index i 
C[i] [s] [t] [n] [k] coefficient C^ n (xk) for the distribution with index i 

S[i] [j] [n] [k] candia.c: coefficient S l n (x k ) for g (j = 0) or (j= 1) 

candia_trunc . c: coefficient S l n {xk) for the distribution with index j 

R[i] [k] value of the PDF f(xk) with index i 



Table 6: Main variables defined in the program. 



const double xtab[] ={le-7, le-6, le-5, le-4 , le-3, le-2,0 . 1 ,0 .3,0 .5,0 .7,0 .9, 1 . }; 
const double Qtab [] ={100 . } ; 

We have also chosen a truncation index k = 6 at NLO and k = 7 at NNLO. For example, the 
command to run the NLO evolution with the truncated method, in this case, is 

candia_trunc .x 1 6 1 1 tl 

while that for the NNLO evolution with the exact method is 
candia.x 2 7 1 1 e2 

We have chosen two representative PDFs: a non-singlet, xq'~\ and a singlet, xg. 

The small difference between the NNLO gluon distributions calculated with the exact and 
the truncated methods are due to the heavy quark matching conditions (I85H87I) . Since a solution 
in closed form is not available for the singlet sector, both the exact and truncated methods solve 
the singlet equations in the same way. By the way, Equations (1851) and f l87l) . whose first-non 
zero contribution is of order a 2 , mix the non-singlet and singlet sectors at NNLO. 

We will therefore concentrate on the differences in the non-singlet sector. At NLO it is 
evident that there is a better agreement for xq^ with Pegasus using the exact method 
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Table 7: Comparison between Candia (exact method) and Pegasus at /zp = hr = 100 GeV 
in the variable flavor number scheme. In each entry, the first number is Candia's result and 
the second one is the difference between Candia's and Pegasus' result. The Les Houches toy 
model at 2 GeV 2 has been used as initial condition. 

(Table [3). At NNLO the truncated method agrees better with Pegasus at low x while the 
exact method agrees better at high x, where the two regions can be approximately separated 
at a value of x in the interval [10 -4 , 10 -3 ]. This peculiar behavior is not surprising, because our 
exact method of solution for the non-singlet at NNLO is based on the closed form solution of 
the DGLAP equation in Mellin-space f l52l . reconstructed via iteration in x-space, as we have 
explained before, and which is not present in PEGASUS. 

In Table [9l where we show the variation of xq(' and xg with the perturbative order, we 
compare the results obtained running Candia with the exact versus the truncated method, 
respectively, having chosen MRST initial conditions. The plots in Figure HJ have been obtained 
evolving the MRST initial conditions from 1 to 100 GeV with Candia (exact method). Notice 
that as we move to NNLO from the LO result, the gluon density, in particular, is drastically 
reduced in the low-x region. For the MRST runs we have modified the file constants. h as 
follows 

#define GRID_PTS 501 
#define NGP 30 
#define ITERATIONS 15 
#define IMTERP.PTS 4 
#define HFT 1 

const double xtab [] ={le-5 , le-4 ,16-3, le-2 ,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.}; 
const double Qtab [] ={100 . } ; 
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Table 8: Same as Table U\ but this time Candia has been run with the truncated method. To 
avoid repetitions, the columns for which the exact and the truncated methods give the same 
results are not shown. 
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Table 9: Comparison between the exact and the truncated method of Candia at /ip = .ur = 
100 GeV in the variable flavor number scheme. In each entry, the first number is the result from 
the exact method and the second one is the difference between the exact and the truncated 
method. The columns for which the exact and the truncated methods give the same results are 
not shown. The MRST parametrizations at 1 GeV have been used as initial conditions. 
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Figure 1: Plots of xq^ and xg at different perturbative orders for Candia (exact method) 
at \ip = fiR= 100 GeV in the variable flavor number scheme. The MRST parametrizations at 
1 GeV have been used as initial conditions. 
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10 Conclusions 



We have presented a new evolution program, Candia, which solves the DGLAP equation with 
high precision and is completely implemented in x-space. We have also briefly discussed the 
types of solutions which are implemented in the program, from the truncated to the exact 
ones. In the non-singlet sector we have shown how to construct the exact solutions analytically 
by the iteration of expressions which resum the simple logarithms of the ratio of the two 
couplings at the initial and final evolution scales. We have also addressed the issue of accuracy 
of the solutions and illustrated the difference between brute force and analytical methods, 
showing the connection between our approach and more traditional approaches based on the 
inversion of Mellin moments. We hope to return in the near future with additional numerical 
implementations, which will provide the user with all the necessary tools so to proceed with an 
independent partonometric analysis of the LHC data on the PDFs. 
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